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Abstract 



CN The positive definite symmetric polymer conformation tensor possesses a unique symmetric square root that satisfies a closed 
evolution equation in the Oldroyd-B and FENE-P models of viscoelastic fluid flow. When expressed in terms of the velocity field 
£3 and the symmetric square root of the conformation tensor, these models' equations of motion formally constitute an evolution 

1—5 in a Hilbert space with a total energy functional that defines a norm. Moreover, this formulation is easily implemented in direct 
numerical simulations resulting in significant practical advantages in terms of both accuracy and stability. 
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1. Introduction 

Familiar models of viscoelastic polymeric fluids present 
challenging problems for both mathematical analysis and nu- 
merical computations. One of the difficulties stems from the na- 
ture of the stress evolution equations. Although there is indeed 
some stress diffusion, the diffusion of polymers is typically or- 
ders of magnitude smaller than for non-polymeric molecules 
and so is often neglected in direct numerical simulations. The 
difficulties manifest themselves both in the form of loss of ac- 
curacy and stability in numerical schemes, and in the absence of 
effective a priori estimates for analysis. Despite recent progress 
in the field, many important problems remain open. 

In this paper we focus on two models, Oldroyd-B and FENE- 
P. It has proven to be a difficult task to devise numerical 
schemes that are efficient, accurate, and stable at the same time. 
One way to ease the numerical problems is to add artificially 
large stress diffusion, and this has been done for a long time. 
Fattal and Kupferman proposed a log-conformation scheme di- 
rectly evolving the matrix logarithm of the positive definite con- 
formation tensor 1 that, according to their reports, indeed helps 
with stability issues. Another method developed by Collins et al 
evolves the eigenvalues of the conformation tensor 2 . Lozinski 
and Owen also proposed to work with the deformation tensor 3 , 
another of the square roots of the conformation tensor. 

We consider a square root method as well, but unlike any 
previous work we are aware of, we derive an evolution equa- 
tion for the positive-definite square root by taking advantage 
of the 0(n) degeneracy in the matrix square root in n dimen- 
sions. This turns out to be both theoretically and numerically 
convenient. On the one hand it allows the dependent variables 



in the Oldroyd-B and FENE-P models to take values in a vec- 
tor space with a natural norm defined by the physical energy. 
On the other hand we observe that, at practically no additional 
computational cost, this formulation produces significant gains 
in both numerical stability and numerical accuracy — without 
adding any artificial stress diffusion — as compared to directly 
evolving the conformation tensor. We note in particular that the 
flow studied in this paper has a hyperbolic stagnation point and 
our simulations for the Oldroyd-B model extend far beyond the 
Weissenberg number at which the stress in the associated steady 
flow becomes infinite. 



2. Mathematical Framework 

The nondimensional equations of motion are 

<9u(x,0 +u>Vu = — Au + V-r+f(x,0, V-u = 0, (1) 
at Re 

with x e K 1 (n = 2 or 3) and Reynolds number Re = Ui/v 
where U and i represent appropriate choices of velocity and 
length scales for the problem under investigation. The exter- 
nally applied body force is denoted by f , and the polymer stress 
tensor r(x, t) is 



Re 



s(c) 



(2) 



where the symmetric positive definite polymer conformation 
(a.k.a. configuration) tensor c(x, t) evolves according to 



<9c T 

— + u • Vc = cVu + (Vu) c + s(c). 

dt 



(3) 
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The parameter s is a coupling constant proportional to the con- 
centration of the polymers in the fluid, and the tensor s(c) 
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takes different forms in various non-Newtonian models. For 
the Oldroyd-B model 



s(c)=-j-(I-c) 
Wi 



(4) 



where the Weissenberg number Wi = UA/£ is the product of 
the polymer relaxation time A and the rate of strain U/€. For 
the FENE-P model 



S(C) " Wif 1 l-(trc// 2 ) 



(5) 



where I 2 is proportional to the square of the maximum poly- 
mer length. For both models the total mechanical energy of the 
system is the sum of the fluid's kinetic energy and the elastic 
potential energy of the polymers: 

S(t) = ^J [|u(x, t)\ 2 + tr t] dx dy dz. (6) 

This energy 1 is formally conserved by the dynamics in the lim- 
its of infinite Reynolds and Weissenberg numbers. 

Unlike the situation for Newtonian fluids modeled by the in- 
compressible Navier-Stokes (or Stokes) equations, the total en- 
ergy does not define a natural norm, or even a metric, in the 
phase space of the dynamical fields u and c. Indeed, the (u, c) 
phase space is not even a linear vector space. This mathemat- 
ical awkwardness results from the fact that the relevant space 
for the conformation tensor c, the space of symmetric positive 
definite matrices, is not a vector space: linear combinations of 
positive matrices are not necessarily positive. These facts com- 
plicate the analysis of these models and preclude implementa- 
tion of useful techniques including nonlinear (energy) stability 
notions 4 . 

This problem can be circumvented, however, by reformulat- 
ing the models in terms of the (unique) symmetric square root 
b(x, t) of the conformation tensor c(x, t). We write 

n 

c ij (x,t) = Y j b ik (x,t)b kj (x,t) with b ij (x,t) = b ji (x,t), (7) 
k=\ 

so the polymer energy density is a function of the matrix norm 
ofb, 



= £fe 2 =tr(b r b) = trc. 



(8) 



i,j=\ 



The work required to implement this proposal is to precisely 
articulate the dynamics of b, a not altogether trivial task due to 
the inherent degeneracy of the matrix square root. 

In the Oldroyd-B case solutions of differential equations of 
the form 



(| + ..v)»- 



bVu + ab+-U(b r )- 1 -b), (9 ) 
zWi 



l 8 does not include the entropic term that contributes to the free energy of 
the system. 



where a(x, t) is any antisymmetric matrix, satisfy b r b = c 
pointwise in space and time when the initial data satisfy 
b r (x,0)b(x,0) = c(x,0). Likewise, in the FENE-P case the 
evolution 



+ u V b 



:bVu+ab+ ^( (bI) 1 



1-1 



II 2 // 2 



, (10) 



produces such a square root of c when a(x, t) is any antisym- 
metric matrix and b r (x, 0)b(x, 0) = c(x, 0). The key observa- 
tion is that by choosing a(x, t) properly we can tune the evolu- 
tion equations (9) and (10) — and similar models with an upper 
convective derivative — to preserve the symmetry of b. That is, 
starting with symmetric initial data b r (x, 0) = b(x, 0), the sub- 
sequent evolution will preserve the symmetry. 
Toward this end we write Vu = (mj) = (up) and 



a = 





-an 
-an 



an 

-a 23 



an 
a 2 3 




, /J =1,2,3, 



(11) 



in n = 3 spatial dimensions and 



012 

-an 



in n = 2 dimensions. Define 



r = (nj) = b(Vu) + ab. 



(12) 



(13) 



We now show that we may choose the matrices a, depending on 
Vu and the symmetric matrix b pointwise in space and time, so 
that r is a field of symmetric matrices, i.e., r t j = rp. For n = 3 
the explicit formulas for the elements come from solving the 
system of 3 linear equations 

(fell + ^22)^12 + ^23^13 - ^31^23 
^23^12 + (^11 + fe33)^13 + fel2^23 
-fel3<2l2 + fel2<3l3 + (fe22 + ^33^23 



Wi, 
W2, 
W 3 , 



(14) 
(15) 
(16) 



where 



W\ 


= (fel2^1,l - 


bnu 2 ,\) 


+ 


(b 22 u h2 


- bnu 2 , 2 ) 








+ 


(fe32^l, 3 


-fe31^2,3)> 


w 2 


= (fel3^1,l - 


fell"3,l) 


+ 


(fe33^1,3 


-fe31^3,3) 








+ 


(fe23^l,2 


-fe21^3,2), 


w 3 


= (fel3^2,l - 


b\ 2 m,\) 


+ 


(fe23^2,2 


- b 22 U3, 2 ) 








+ 


(fe33^2,3 


- fe32^3,3). 



In matrix notation, this is the system of equations 



( fell +fe22 fe23 

fe23 fell+fe33 
-fe 3 l fei2 



-fe 3 l 
fel2 
fe22 + fe3 





' a l2 * 








an 




w 2 


/ 


< a 23 , 







(17) 
(18) 
(19) 

(20) 



Then by swapping the first and the third columns of the coef- 
ficient matrix (and hence, also a 2 3 and a 12), and subsequently 
swapping the first and the third rows of the resulting coefficient 
matrix (and hence, also w\ and w 3 ), and finally multiplying the 



2 



second row and the second column of the resulting matrix by 
-1 (and hence, also replacing and w 2 by -an and -w 2 , re- 
spectively), we obtain 



where 





(tr(b)I 


-b) a = 


= v, 




' a 23 " 




W3 


a = 




v = 


-W2 




< a l2 , 







(21) 



(22) 



When b is symmetric at the point (x, t) there is an orthogonal 
matrix p(x, t) such that 

b = p r diagUi,^ 2 ,^3}p, (23) 
where X[ are the eigenvalues of b. Thus, we have 

trace(b)I - b = trace(b)I - p r diag {X\ , A 2 , h} P (24) 
= p r (trace(b)I - diag {A u A 2 , ^ 3 })p (25) 
= p r diag {/l 2 + /I3, /ii + A3, X\ + /l 2 } p (26) 

Again, assuming that b is positive definite (although this condi- 
tion can be clearly relaxed to include a large class of semidef- 
inite objects) we can solve for a uniquely so that the evolution 
equations (9) and (10) used to obtain b at later times are sym- 
metrized. The explicit algebraic formulas for the elements 
for n = 3 are displayed in the appendix. In the much simpler 
case ofn = 2 space dimensions we have 



an 



(bnu\,\ ~ b n u 2 ,i) + (b 22 u h2 - b 2 iu 2 , 2 ) 



fen +b 



(27) 



'22 



This construction puts the full dynamics back in a vector 
space setting (the direct product of vector fields u and sym- 
metric tensor fields b). For Oldroyd-B the energy functional 
(6) is proportional to the vector norm (squared) on the direct 
product space (modulo an additive constant). For the FENE-P 
model, the convexity of b r b/(l - (tr(b r b)// 2 )) in a neighbor- 
hood of the origin (i.e., for ||b|| < /) allows as well for a natural 
energy norm. And as shown in the next section, in some cases 
this reformulation of the dynamics leads to significant practical 
improvements in direct numerical simulations. 

3. Numerical Experiments 

As a test of the numerical accuracy and stability of the 
square-root method we consider the zero-Reynolds number 
(Stokes) limit of the Oldroyd-B and FENE-P models in which 
case the momentum equation (1) reduces to 



Vp = Au + V • r + f , V • u = 0. 



(28) 



Here r = -ss(c) with s(c) given by (4) or (5) for the Oldroyd-B 
model or FENE-P model respectively. In what follows we fix 
s = 0.5. Following recent studies 5 ' 6 we consider a 2^-periodic 
domain in n - 2 space dimensions ([-7r,7r] 2 ) and impose a 
steady background force 




f = (-2 sin x cos y, 2 cos x sin y), 



(29) 



Figure 1: (a) Contour plots of curlf for the force given by Eq. (29). (b) - (d) 
Contour plot of vorticity, trc, and C12 for isotropic initial data, Wi = 5 at t = 10. 



curlf is shown in Fig. 1(a). In the absence of polymer stress this 
yields a four-roll mill geometry for the fluid velocity. One well- 
known consequence of this body-force imposed extensional ge- 
ometry in the Oldroyd-B model is that the polymer stress and 
stress gradients grow exponentially in time 5,7,8 and inevitably 
produce numerical problems. In particular, when resolving 
steep gradients the loss of positive definiteness of the confor- 
mation tensor due to numerical error can lead to breakdown 
of the computational schemes. One common solution to these 
difficulties is the addition of artificial polymer stress diffusion. 
Although some polymer stress diffusion can be derived from 
the basic physics in the model, the magnitude of the physically 
relevant diffusion is far too small to have an effect on numer- 
ical simulations 9 . In the following we do not add any stress 
diffusion to the numerical calculations. 

The Stokes-Oldroyd-B system (and FENE-P) is solved with 
a pseudo-spectral method 10 . In the usual formulation the 
conformation tensor is evolved using a second-order Adams- 
Bashforth method. The initial data is prescribed, and given c 
the Stokes equation is inverted in Fourier space for u. Given 
u the nonlinearities of the conformation tensor evolution are 
evaluated using a smooth filter applied in Fourier space before 
the quadratic terms are multiplied in real space 11 . The con- 
formation tensor is then discretized on the Fourier transform 
side. It should be noted that the numerical implementation of 
the square-root method adds no computational cost. 

In a recent investigation 5 the Stokes-Oldroyd-B equations 
were solved starting from isotropic, i.e., c(0) = I, initial data 
and the stress was observed to diverge (exponentially in time) 
at the extensional stagnation points in the flow for sufficiently 
large Weissenberg number. However, outside of an exponen- 
tially decreasing region around the extensional stagnation point, 
the solutions became steady after an initial transient. These 
near- steady solutions preserve many symmetries: the stress is 
symmetric and aligned along the direction of extension and the 
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flow has an underlying four-roll structure. For sufficiently large 
Wi additional oppositely signed vortices arise along the stable 
and unstable manifolds of the extensional stagnation point. Fig- 
ure 1(b) - (d) displays these symmetric solutions in the case 
Wi = 5, at t = 10. These symmetries are broken as the initial 
data is perturbed and it was shown that instabilities arise for 
sufficiently large Wi 6 ' 12 . 

In what follows we discuss both accuracy and stability im- 
provements for the Oldroyd-B and FENE-P models using the 
square root method. In subsection 3.1 we consider homoge- 
neous isotropic initial data c(x, 0) = b 2 (x, 0) = I and compare 
solutions to the Stokes-Oldroyd-B model obtained using the di- 
rect evolution of c with those obtained by evolving the symmet- 
ric square root. In subsection 3.2 we consider perturbations to 
the initial data for the Stokes-Oldroyd-B model (as in previous 
studies 6,12 ) to see how far the simulations run, for a fixed reso- 
lution using each method, before numerical divergences appear. 
Finally, in subsection 3.3 we revisit both of these questions for 
the FENE-P model. 

3.1. Accuracy 
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Figure 2: (a) Absolute error |c^ -c exa ctl measured along the axis of compression 
of the first component of the conformation tensor with N 2 = 256 2 , Wi = 5, at 
t = 10. Dotted line is c and solid line is b 2 . (b) Relative error |Cj y~ Cex f 1 , 

V J ICexactl 





Figure 3: (a) Relative error in L l shows 
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computed with c and 



b , for Wi = 5 at t = 10. (b) Improvement in accuracy as a function of N, for 
Wi = 1 - 5 at T = f/Wi = 2. 



Figure 3 (a) shows the relative error 

IjCjy - CexactllL 1 
IICexactllL 1 



(30) 



measured in the L l ([-n, n] 2 ) norm for both and (b 2 )^v for 
N 2 = 32 2 ,64 2 , 128 2 ,256 2 ,512 2 , comparing to the "exact" so- 
lution defined by the N 2 = 2048 2 simulation. The L 1 norm is 
chosen because it takes into account the average error over the 
entire domain. This computation is also for Wi = 5 at t = 10. 

In this averaged sense we see that the error is always smaller 
using the square root method. The improvement in accuracy (in 
the L 1 -sense) using the square root method is shown in figure 
3 (b). Here we plot 



|error c - error b 2 1 
|error c | 



(31) 



forWi = 1, 2, 3,4, 5 for N 2 = 32 2 , 64 2 , 128 2 , 256 2 , 512 2 . Each 
simulation is computed at T = t/Wi = 2, and this scaled-time is 
used because the solutions grow exponentially like e t/Wl . There 
is a significant improvement for higher Wi, in particular for 
lower resolutions N 2 = 128 2 and 256 2 . 



Figure 2 shows the difference between the solution to Stokes- 
Oldroyd-B with N 2 = 25 6 2 grid cells and the "exact" solution 
computed by evolving c with N 2 = 2048 2 grid cells, resolved to 
at least 6 digits of accuracy. The dotted line shows the solution 
computed by directly evolving the conformation tensor c, while 
the solid line shows the solution computed by evolving the sym- 
metric square root. The simulation is performed with Wi = 5 
and the result of the computation is shown at t = 10. Panel (a) 
shows the absolute error (|c# - Cexactl) m me ^ rst com P onent 
of the conformation tensor (en) along the direction of com- 
pression because this is precisely where steep gradients form. 
We observe that away from the extensional stagnation point the 
square root gives a much better approximation and it is only 
very near the extensional stagnation point where the evolution 
of c gives a better approximation. The relative error is shown in 
panel (b) to emphasize that although the square-root method's 
error at the extensional stagnation point is larger, the relative er- 
ror is actually quite small because the stress is very large there. 



3.2. Stability 

Experiments on low-Reynolds number viscoelastic turbu- 
lence 13 ' 14 ' 15 and instabilities in extensional geometries 16 have 
inspired many numerical studies of low-Reynolds number vis- 
coelastic fluids 17 ' 18 ' 19 ' 6 ' 12 . Two main instabilities were observed 
in one investigation 6,12 : first, for sufficiently large Wi, if a small 
perturbation is introduced in the initial stress conformation the 
extensional stagnation point in the flow becomes unstable and 
loses the pinning to the background steady force. For larger 
Wi other stagnation points lose their pinning to the background 
force and higher oscillations arise in the flow. These instabili- 
ties occur on long time scales and some artificial polymer stress 
diffusion was introduced to fully resolve the stress 6 ' 12 . Here we 
test these same perturbations to the initial data and for a fixed 
resolution N 2 = 25 6 2 and run the same simulations without any 
artificial polymer stress diffusion. 

Figure 4 (a) shows a plot of the first component of the confor- 
mation tensor for Wi = 10 at t = 15 computed both by evolving 
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Figure 4: (a) Plot of first component of conformation tensor cn(0,y) and 
b^ 1 (0,_y) along direction of compression near the extensional stagnation point 
for Wi = 10 at t = 15. Computations of c stop producing finite values at t = 19. 
(b) Plot of max(tr(b 2 )) as a function of time over < t < 1500. (c) Wi = 10, 
contour plot of tr(b 2 ) at t = 1000 on [-7r,7r] 2 . (d) Wi = 10, contour plot of the 
vorticity of the flow field at t = 1000 on [-n, -n] 2 . All simulations performed 
with N 2 = 256 2 grid points. 



the conformation tensor (solid line) and the square root (dotted 
line). The plot is shown along the axis of compression and it is 
evident that the stress has accumulated significantly and is quite 
large near the extensional stagnation point in the flow (y = 0). 
The oscillations produced in the direct evolution of c lead to the 
loss of positive-definiteness of the conformation tensor, and the 
numerical scheme breaks down. This figure shown is at t = 15 
and the computation fails to produce finite numbers at t = 20. 
However using the square root one can run these simulations 
to t = 1500 and even beyond. Figure 4 (b) shows a plot of 
max(tr(b 2 )) as a function of time for < t < 1500. It is impor- 
tant to point out that although max(tr(b 2 )) remains bounded in 
this case (with no artifical stress diffusion) this quantity clearly 
depends on N and this is one way that the accuracy of the so- 
lution is lost. However, it is not clear that this level of loss 
of accuracy is entirely relevant to the flow because the region 
where tr(b 2 ) gets large diminishes exponentially in time even as 
tr(b 2 ) grows exponentially in time 5 . Figure 4 (c) and (d) show 
contour plots of tr(b 2 ) and the vorticity of the flow on [-n, n] 2 . 
Here we see that the previously observed instabilities 6 ' 12 are at 
least qualitatively reproduced. The time-dependent behavior is 
similar, too: the four-roll mill structure of the background force 
is preserved initially, the extensional stagnation point leaves the 
origin, and eventually time-dependent oscillations arise in the 
flow. 

Of course with fixed resolution and no stress diffusion there 
is an inevitable loss of accuracy. This can be seen in Fig. 4 
(c) and (d) in the slightly fuzzy images indicating oscillations 
while attempting to resolve the steep gradients in the conforma- 
tion tensor and vorticity. It is noteworthy that these simulations 
are performed with no artificial stress diffusion but nevertheless 



qualitatively reproduce the well-resolved results that utilized ar- 
tificial diffusion 6,12 . The same computations simply cannot be 
performed with a direct evolution of the conformation tensor 
(in this particular code). The square root method allows sim- 
ulations to run much longer and at much higher Weissenberg 
number than evolving c directly allows. This indicates that it 
might be possible to use much smaller — closer to the physi- 
cally realistic quantity — stress diffusion and still obtain reason- 
ably accurate results, although this will not be pursued in this 
paper. 

3.3. FENE-P 



Relative error 
along axis of compression 



Relative error in L 1 




Figure 5: (a) Relative error measured along the axis of compression 

of the first component of S(c) for FENE-P with N 2 = 256 2 , Wi = 5 at t = 10, 
I 2 = 100. Dotted line is S(c) and solid line is S(b 2 ). (b) Relative error in L 1 

shows l|S ^ ~ Se ^ ctllLl computed with S(c) and S(b 2 ), for Wi = 5 at t = 10, 

I Inexact | \ L 1 

I 2 = 100. (c) Wi = 20, contour plot of tr(S(b 2 )) at t = 100 on [-n,n] 2 , I 2 = 
225. (d) Wi = 20, contour plot of the vorticity of the flow field at t = 100 on 
[-n, -n] 2 , I 2 = 225. All simulations done with N 2 = 256 2 grid points, (e) 
Wi = 50, contour plot of tr(S(b 2 )) at t = 500 on [-tt, tt] 2 , I 2 = 225. (f) Wi = 50, 
contour plot of the vorticity of the flow field at t = 500 on [-n, -n] 2 , 1 2 = 225. 
All simulations done with N 2 = 25 6 2 grid points. 



The FENE-P model enforces a limit (I 2 ) on the magnitude of 
tr c so the conformation tensor remains bounded. Steep gradi- 
ents can still arise in the polymer stress, however, and numeri- 
cal difficulties remain. Therefore we also simulated the FENE- 
P model in a Stokesian solvent to check for possible accuracy 
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and stability improvements by evolving the symmetric square 
root. 

Figure 5 (a) and (b) are analogs of Figures 2 (b) and 3 (a) for 
FENE-P. The simulations were performed with Wi = 5 and cut- 
off I 2 = 100, and are displayed at t = 10. Rather than plot the 
conformation tensor c and b 2 , however, it is more analogous to 
plot S = 1-(t rc //2) because this is closely related to the physical 
stress tensor and includes the factor that gets very large as trc 
gets near the cut-off I 2 . The accuracy improvement is not as 
large here as it was for Oldroyd-B: for Wi = 5 the improvement 
is about 65% at N = 256 2 and is only 31% for N 2 = 512 2 
but there is still some improvement (especially away from the 
extensional stagnation point). As before, the "exact" solution 
here comes from a simulation with N 2 = 1024 2 . 

The significance of the symmetric square root method for 
FENE-P is much more apparent in terms of stability. The fact 
is that we can increase Wi much more than we can for Stokes- 
Oldroyd-B. We show results from two different simulations to 
demonstrate this. First in Fig. 5 (c) and (d) we show results 
from perturbed initial data with Wi = 20 at t = 100, with 
I 2 = 225. This is just after the onset of the instability and the 
flow is still nearly symmetric. The stress has accumulated along 
the incoming and outgoing streamlines of the extensional stag- 
nation point and the four-roll mill structure of the vorticity is 
still largely preserved. Fig. 5 (e) and (f) show results from per- 
turbed initial data with Wi = 50 at t = 500, with I 2 = 225. 
These same computations evolving c fail to produce finite val- 
ues before t = 60 whereas the evolution of b appears to con- 
tinue indefinitely — although, again, there must be some loss of 
accuracy. The qualitative behavior is similar to the solutions of 
Stokes-Oldroyd-B and the instabilities discussed for that case 
also occur here. In Fig. 5 (c) and (e) we show contour plots of 
tr(S(b 2 )) after the instabilities have developed and observe that 
max tr(S) is quite large. The time-dependent behavior is also 
quite complicated and as one can see in Fig. 5 (f), the vorticity 
of the flow is also very complex with many additional vortices 
continually arising and being destroyed in the flow. 

4. Discussion and Conclusions 

In hindsight both the symmetrization procedure and directly 
computing the square root evolution equations (9) and (10) 
might have been expected to contribute to the gains in stability 
and accuracy. Taking the square root reduces large amplitudes 
which, not unexpectedly, reduces the stiffness in time stepping. 2 
Moreover symmetrizing the system may reduce the stiffness of 
the time marching as compared to taking a = and simply 
computing the deformation tensor because components of the 
symmetric square root matrix will generally have less variance 



than a square root with no symmetry. And computing the square 
root instead of c has other advantages: the square root computa- 
tion ensures positivity of c in the numerical scheme compared 
to the most direct evolution of the conformation tensor. We 
have observed that in practice the square root method can be 
applied at higher Wi and for longer time without any artificial 
numerical stress diffusion than evolving the conformation ten- 
sor directly can, enabling one to obtain numerical solutions in 
more situations. 

One less obvious but perhaps important advantage is the fol- 
lowing. Assuming the locality of modal interactions, i.e., that 
lower spectral modes of c are determined predominantly by 
lower modes of the square root matrix b, we can expect good 
information about up to the first 2N modes of the conformation 
tensor c when we know just the first N modes of the square root 
b. This speculation basically boils down to the assumption that 
the Galerkin approximation method works well for both b and 
c for sufficiently large N. Thus we might expect that evolution 
of the square root improves both stability and accuracy, at least 
in spectral or pseudo spectral schemes. Whether this is the case 
in other types of spatial discretizations requires further investi- 
gation. Of course it also remains to implement the full n = 3 di- 
mensional symmetric square root algorithm and systematically 
compare its performance with conventional schemes used to in- 
vestigate, for example, turbulent drag reduction 20 . 

We emphasize that all the numerical simulations shown have 
been performed without artificial diffusion. One can add artifi- 
cial diffusion to the advection of b, however for this to match 
the standard diffusion of the conformation tensor, a more com- 
plicated nonlinear diffusion term is needed. The square-root 
method still has limitations for sufficiently large Wi so it is nat- 
ural to ask if one could use a more physically realistic diffusion 
coefficient with the square-root method and this will be pursued 
in future work. 

The advantage of expressing the Oldroyd-B and FENE-P 
models as evolutions in a vector space (indeed, a Hilbert space) 
remains one of theoretical elegance at this point. Whether or not 
this formulation might assist the rigorous mathematical analysis 
of these models is an open question. 
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Appendix A. Entries of the antisymmetric matrix for n — 
3 

The entries a \ 2 , a 13, a 23 of the antisymmetric matrix a in (1 1) 
are given by 

Da n = (TiT 2 -Bl)wi - {B X T X + B 3 B 2 ) w 2 

+ (B 2 T 2 + B l B 3 )w 3 , (A.l) 

Da l3 = -(BiTi+B 3 B 2 )wi + (t x T 3 - B 2 2 )w 2 

- (B 2 B x +B 3 T 3 )w 3 , (A.2) 

Da 23 = (B 2 T 2 + BiB 3 )wi - (B 2 B X + B 3 T 3 )w 2 

+ (r 2 r 3 -^)w 3 , (A.3) 

where 

D = det((trb)I-b) 

= T x (T 2 T 3 - Bj) - B 2 (B 2 T 2 + B X B 3 ) 

- B 3 (B 2 B l +B 3 T 3 ), (A.4) 

Ti = b 22 + b 33 , T 2 = b n + b 33 , T 3 = b n + ^22, (A.5) 

B 3 = &12, #2 = bis, Bi = £23, (A.6) 

and w\,w 2 , and W3 are given in (17)-(19). 
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